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ABSTRACT 

When  one  considers  the  eflFect  in  the  physical  space,  Daubechies-based  wavelet 
methods  are  equivalent  to  finite  difference  methods  with  grid  refinement  in  regions  of 
the  domain  where  small  scale  structure  exists.  Adding  a  wavelet  basis  function  at  a 
given  scale  and  location  where  one  has  a  correspondingly  large  wavelet  coefficient  is, 
essentially,  equivalent  to  adding  a  grid  point,  or  two,  at  the  same  location  and  at  a 
grid  density  which  corresponds  to  the  wavelet  scale.  This  paper  introduces  a  wavelet- 
optimized  finite  difference  method  which  is  equivalent  to  a  wavelet  method  in  its 
multiresolution  approach  but  which  does  not  suffer  from  difficulties  with  nonlinear 
terms  and  boundary  conditions,  since  all  calculations  are  done  in  the  physical  space. 
With  this  method  one  can  obtain  an  arbitrarily  good  approximation  to  a  conservative 
difference  method  for  solving  nonlinear  conservation  laws. 
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1  Introduction 


In  the  numerical  simulation  of  equations  which  model  physics  it  is  common  that  small 
scale  structure  will  exist  in  only  a  small  percentage  of  the  domain.  If  one  chooses 
a  uniform  numerical  grid  fine  enough  to  resolve  the  small  scale  features  then  in  the 
majority  of  the  domain  the  solution  to  the  equations  will  be  over  resolved.  One  would 
like,  ideally,  to  have  a  dense  grid  where  small  scale  structure  exists  and  a  sparse  grid 
where  the  solution  is  composed  only  of  large  scale  features. 

Consider  now  a  Daubechies-based  wavelet  system.  Wavelets  provide  a  natural 
mechanism  for  decomposing  a  solution  into  a  set  of  coefficients  which  depdhd  on  scale 
and  location.  One  can  then  work  with  the  solution  in  a  compressed  form  where  one 
works  only  with  the  wavelet  coefficients  which  are  larger  in  magnitude  than  a  given 
threshold.  Wavelets,  therefore,  soimd  ideal  for  solving  the  type  of  problem  mentioned 
the  previous  paragraph.  There  are,  however,  serious  problems  matching  the  order 
of  differentiation  accuracy  at  the  boundary  for  nonperiodic  boundary  conditions,  see 
[9],  with  the  superconvergence  encountered  with  periodic  boundary  conditions,  see 
[7].  Furthermore,  wavelet  methods  generally  require  a  tranformation  between  the 
physical  space  and  the  coefficient  space  for  either  evaluation  of  nonlinear  terms  or  for 
differentiation. 

In  this  paper  a  wavelet  method  which  satisfies  the  goals  of  the  first  paragraph 
while  using  the  wavelet  machinery  outlined  in  the  second  paragraph  without  the 
accompanying  difficulties  encountered  at  the  boundaries  and  the  expense  of  constantly 
tranforming  between  the  physicsd  space  and  the  coefficient  space  will  be  introduced. 
That  is,  the  new  method  utilizes  the  strength  of  wavelets,  scale  detection  and  data 
compression,  while  avoiding  the  difficulties  by  using  wavelets  in  their  finite  difference 
form. 

The  following  is  a  list  of  the  sections  of  this  paper  with  the  noteworthy  points. 
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1.  Introduction 


2.  Wavelet  Definitions  and  Notation 

3.  Finite  Difference  Grid  Refinement  and  Wavelets: 

This  section  will  establish  that  Daubechies-based  wavelet  methods  are  equal  to 
finite-difference  methods  with  grid  refinement. 

4.  The  Wavelet-Optimized,  Adaptive  Grid,  Finite  Difference  Method: 

In  this  section  a  new  numerical  method  which  utilizes  the  strength  of  wavelets 
and  avoids  the  difBculties  will  be  proposed.  That  is,  wavelets  will  be  utilized 
to  define  the  grid  for  finite  difference  methods.  The  new  method  is  named, 
*The  Wavelet-Optimized,  Adaptive  Grid,  Finite  Difference  Method’,  or  simply 
‘WOFD’. 

5.  WOFD  applied  to  Burgers*  equation: 

Numerical  results  of  the  wavelet-optimized,  adaptive  grid,  finite  difference  method 
applied  to  Burgers’  equation  with  periodic  and  nonperiodic  boundary  conditions 
will  be  given. 

6.  Accuracy  of  WOFD: 

The  error  in  finite-difference  derivative  approximation  on  a  5-point  stencil  is  of 
the  form. 

Err  =  A,AjA3A4i/W(o). 

Think  of  /(x)  as  a  pure  mode,  /(x)  =  e’*^,  where  (  is  frequency  or  wave  number. 
When  the  data  is  locally  smooth,  i.e.,  composed  of  low  frequencies,  the  wavelet 
coefficients  are  small  and  consequently  the  A’s  are  allowed  to  be  large.  When 
the  data  is  locally  oscillatory,  i.e.,  composed  of  high  frequencies,  the  wavelet 
coefficients  are  large  and  WOFD  reduces  the  size  of  the  A’s.  The  effect  is 
that  the  derivative  approximation  error  will  not  grow  faster  than  linearly  with 
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respect  to  frequency.  Recall  that  without  grid  adjustment  this  error  would  grow 
as  a  fifth  power  of  frequency  for  a  fourth  order  scheme. 

7.  Stability  of  WOFD: 

Analytical  stability  methods  are  beyond  reach  due  the  arbitrary  nature  of  the 
grid.  But,  in  practice  the  method  displays  no  instability  when  applied  to  Burg¬ 
ers’  equation. 

8.  Efficiency  of  WOFD: 

The  WOFD  method  finds  an  approximation  to  the  solution  found  on  the  finest 
scale  across  the  whole  domain.  The  efficiency  depends  on  the  rate  of  data 
compression.  That  is,  if  the  finest  scale  has  N  grid  points  and  the  WOFD 
averages,  say.  No  grid  points,  then  the  WOFD  method  will  find  the  solution 
using,  roughly,  ^  times  the  number  of  operations  used  to  find  the  finest  grid 
solution. 

9.  WOFD  in  Higher  Dimensions 

The  discussion  here  will  be  limited  to  grid  selection.  It  will  be  seen  from  a  few 
examples  that  WOFD  is  an  effective  method  for  grid  selection  in  higher  dimen¬ 
sions.  The  examples  given  are  for  two  dimensions.  The  local  ‘spectral  analysis’ 
of  a  wavelet  method  provides  exactly  the  information  needed  to  thoroughly 
understand  the  data  and,  hence,  define  a  grid  properly. 

10.  Conclusion: 

The  WOFD  method  is  an  efficient  and  stable  alternative  to  a  Daubechies  wavelet 
method.  The  WOFD  method  and  a  wavelet  method  are  essentially  the  same. 
The  only  significant  difference  is  the  manner  in  which  the  grid  is  refined.  The 
WOFD  method,  by  contrast,  avoids  difficulties  with  nonlinear  terms  and  bound¬ 
aries  by  performing  all  calculations  in  the  physical  space. 


3 


2  Wavelet  Definitions  and  Relations 


The  term  wavelet  is  used  to  describe  a  spatially  localized  function.  ‘Localized’  means 
that  the  wavelet  has  compact  support  or  that  the  wavelet  almost  h2is  compact  sup¬ 
port  in  the  sense  that  outside  of  some  interval  the  amplitude  of  the  wavelet  decays 
exponentially.  We  will  consider  only  wavelets  that  have  compact  support  and  that 
are  of  the  type  defined  by  Daubechies  [4]  which  are  supported  on  [0,2 A/  —  1],  where 
M  is  the  number  of  vanishing  moments  defined  later  in  this  section. 

To  define  Daubechies  wavelets,  consider  the  two  functions  (f>{x)  and  tj){x)  which 
are  solutions  to  the  following  equations: 

<l,{x)  =  V2^^hk<f>{2x-k),  (1) 

ib=0 

L-l 

ip{x)  =  >72  gk<l>{2x  -  k),  (2) 

fc=0 

where  <f){x)  is  normalized, 

r  <i>{x)dx = 1.  (3) 

J — oo 

Let, 

=  2-i^l2-’x  -  k),  (4) 

and 

«(!)  =  2-^(2-‘x  -  k),  (5) 

where  j,  k  €  Z,  denote  the  dilations  and  translations  of  the  scaling  function  and  the 
wavelet. 

The  coefficients  H  =  and  G  =  are  related  by  gk  =  (— l)^AL_jt  for 

A:  =  0, ...,  L  —  1.  Furthermore,  H  and  G  are  chosen  so  that  dilations  and  translations 
of  the  wavelet,  form  an  orthonormal  basis  of  L^{R)  and  so  that  V’(®)  has  M 

vanishing  moments.  In  other  words,  ^^(x)  will  satisfy 


where  Su  is  the  Kronecker  delta  function.  Also,  ij}{x)  =  V’o(^)  satisfies 

f  0(x)x”‘dar  =  0,  (7) 

J—OO 

for  m  =  0,  —  1.  Under  the  conditions  of  the  previous  two  equations,  for  any 

function  /(x)  €  L>^{R)  there  exists  a  set  {djt}  such  that 

/(*)  =  (8) 

j€ZKeZ 

where 

djk=  [  f{x)tf)i{x)dx.  (9) 

J— oo 

The  number  of  vanishing  moments  of  the  wavelet  V’(®)  defines  the  accuracy  of 
approximation.  The  two  sets  of  coefficients  H  and  G  are  known  in  signal  processing 
literature  as  quadrature  mirror  filters  [5].  For  Daubechies  wavelets  the  number  of 
coefficients  in  H  and  G,  or  the  length  of  the  filters  H  and  G,  denoted  by  L,  is  related 
to  the  number  of  vanishing  moments  M  by  2M  =  L.  For  example,  the  famous  Haar 
wavelet  is  found  by  defining  H  as  ho  =  hi  =  1.  For  this  filter,  H,  the  solution  to 
the  dilation  equation  (1),  ^(x),  is  the  box  function:  ^(x)  =  1  for  x  €  [0,1]  and 
^(x)  =  0  otherwise.  The  Haar  function  is  very  useful  as  a  learning  tool,  but  it 
is  not  very  useful  as  a  basis  function  on  which  to  expand  another  function  for  the 
important  reason  that  it  is  not  differentiable.  The  coefficients,  H,  needed  to  define 
compactly  supported  wavelets  with  a  higher  degree  of  regularity  can  be  found  in  [4]. 
As  is  expected,  the  regularity  increases  with  the  support  of  the  wavelet.  The  usual 
notation  to  denote  a  Daubechies  wavelet  defined  by  coefficients  H  of  length  L  is  Di¬ 
li  is  usual  to  let  the  spaces  spanned  by  ^(x)  and  ^^(x)  over  the  parameter  k, 
with  j  fixed,  to  be  denoted  by  Vj  and  Wj  respectively: 


The  spaces  Vj  and  Wj  are  related  by  [4] 


...  C  Vi  C  Vo  C  V_,  C  ..., 


(12) 


and  that 


Vi  =  Vrt.0»V+..  (13) 

The  previously  stated  condition  that  the  wavelets  form  an  orthononnal  basis  of  L^{R) 
can  now  be  written  as, 

L\R)  =  ®Wi.  (14) 

iez 

Two  final  properties  of  the  spaces  Vj  are  that 


n  ''i = {Oh  (15) 

J€2 


and 


U  =  l-'W-  (16) 

ieZ 

Of  cotirse,  infinite  sums  and  unions  are  meaningless  when  one  begins  to  implement 
a  wavelet  expansion  on  a  computer.  In  some  way  one  must  limit  the  range  of  the 
scale  parameter  j  and  the  location  parameter  k.  Consider  first  the  scale  parameter  j. 
As  stated  above,  the  wavelet  expansion  is  complete:  L^{R)  =  Therefore, 

any  f{x)  €  L^iR)  can  be  written  as. 


/(®)  =  E  X) 

jezk^z 

where  due  to  orthonormality  of  the  wavelets  /(x)V’i(a))-  In  this  expan¬ 

sion,  functions  with  arbitrarily  small-scale  structures  can  be  represented.  In  practice, 
however,  there  is  a  limit  to  how  small  the  smallest  structure  can  be.  This  would 
depend,  for  example,  on  how  fine  the  grid  is  in  a  numerical  computation  scenario  or 
perhaps  what  the  sampling  frequency  is  in  a  signal  processing  scenario.  Therefore, 
on  a  computer  an  expansion  would  take  place  in  a  space  such  as 


Vo  =  IVi  ©  Wa  0  ...  ©  WO  e  V;, 


(17) 
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and  would  appear  as, 


kez  i=l k€Z 

where  again  due  to  orthonormality  of  the  basis  functions  /(a:)V’fe(a:),  and 

s{  =  f{x)<l>'l{x).  In  this  expansion,  scale  j  =  0  is  arbitrarily  chosen  as  the  finest 

scale  that  is  needed,  and  scale  J  would  be  the  scale  at  which  a  kind  of  local  average, 
(f>i{x),  provides  sufficient  large  scale  information.  In  language  that  is  likely  to  appeal 
to  the  electrical  engineer  it  can  be  said  that  4i{x)  represents  the  direct  current  portion 
of  a  signal  at  location  k  and  that  tl’Kx)  represents  the  alternating  current  portion  of 
a  signal  at,  very  roughly,  frequency  j  and  location  k.  As  stated  above,  one  must 
also  limit  the  range  of  the  location  parameter  k  If  one  assumes  periodicity,  then  the 
periodicity  of  f{x)  induces  periodicity  on  all  wavelet  coefficients,  and  d].,  with 
respect  to  k.  Without  periodicity,  the  location  parameter  k  will  begin  at  1  with  the 
left-hand  side  boundary  functions  and  end  with  some  maximum  number  N  at  the 
right-hand  side  boundary  functions. 

This  completes  the  definition  of  wavelets. 
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3  Finite  Difference  Grid  Refinement  and  Wavelets 


In  this  section  it  will  be  shown  that  Daubechies-based  wavelet  methods  when  con¬ 
sidered  in  the  physical  space  are  equivalent  to  explicit  finite  difference  methods  with 
grid  refinement.  In  a  Daubechies  wavelet  method  the  ‘refinement’  is  accomplished  by 
adding  wavelet  bases  functions  in  regions  where  structure  exists  corresponding  to  the 
scale  of  the  wavelet  used  for  analysis.  In  a  finite  difference  method  the  ‘refinement’  is 
accomplished  by  adding  grid  points  in  regions  chosen  by  some  grid  refinement  mech¬ 
anism.  In  this  section  it  is  argued  that  since  wavelet  methods  correspond  to  central 
finite  difference  operators  when  the  grid  is  uniform  and  since  wavelet  methods  contain 
a  natural  and  effortless  mechanism  for  ‘grid  refinement’,  then  one  can  simply  use  the 
wavelets  to  refine  a  grid  for  finite  difference  operators.  In  this  way  one  can  maintain 
the  superconvergence  encountered  with  periodic  boundary  conditions,  see  [7],  which 
is  lost  when  one  constructs  wavelets  on  an  interval,  see  [9].  That  is,  boundary  condi¬ 
tions  are  imposed  in  the  same  manner  as  for  finite  difference  operators.  Furthermore, 
there  is  no  longer  a  diflSculty  with  nonlinear  terms  requiring  constant  transformation 
between  the  physical  space  and  the  coefficient  space  since  all  calculations  are  done  in 
the  physical  space. 

This  section  contains  four  subsections: 

1.  The  wavelet  decomposition  matrix  will  be  constructed. 

2.  It  will  be  seen  that  under  the  assumption  of  periodicity  and  without  data  com¬ 
pression  that  the  effect  in  the  physical  space  of  differentiation  in  the  D4  subspace 
Vo  is  exactly  the  same  as  differentiation  with  the  optimal  central  4th-order  finite 
difference  operator. 

3.  Now  we  compare  wavelets  and  finite  difference  in  the  subspace  Vo  =  Wi  0  Vi .  If 
Ax  is  the  grid  spacing  in  Vo  then  2Ax  is  the  grid  spacing  in  Vi  and  the  wavelet 
coefficients  in  W]  indicate  if  refinement  is  needed  for  a  local  grid  spacing  of  Ax. 
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4.  Finally,  a  division  of  the  subspace  Vq  which  might  be  used  in  practice  is  studied: 
Vo  =  Wi  ®  Wj  ©  W3  0  V3.  Similar  to  above,  the  grid  spacing  in  the  subspace  V3 
would  be  8Ax.  The  first  refinement  is  controlled  by  the  subspace  W3  which  can 
refine  locally  to  a  grid  spacing  4Aa;.  Likewise,  the  subspace  W-i  refines  locally 
to  a  grid  spacing  of  2Ax  and  Wx  to  a  local  grid  spacing  of  Az. 


3.1  Wavelet  Decomposition  Matrix 

The  wavelet  decomposition  matrix  is  the  matrix  embodiment  of  the  dilation  equation 
defining  the  scaling  function  and  the  accompanying  equation  defining  the  wavelet. 
The  following  two  recursion  relations  for  the  coefficients  and  d\  can  be  found  from 
equations  (1)  and  (2),  respectively: 

n=2M 

4=  E  (19) 

n=l 


~  9n^n+2k-2'  (20) 

n=l 

Denote  the  decomposition  matrix  embodied  by  these  two  equations,  assuming  peri¬ 
odicity,  by  where  the  matrix  subscripts  denote  the  size  of  the  matrix,  and  the 

superscripts  indicate  that  P  is  decomposing  from  scaling  function  coefficients  at  scale 
j  to  scaling  function  and  wavelet  function  coefficients  at  scale  j  -1- 1.  Let  Sj  contain 
the  scaling  function  coefficients  at  scale  j.  (Note  that  when  vector  notation  is  used 
the  scale  is  given  as  a  subscript.)  P  therefore  maps  Sj  onto  Sj+i  and  dj+i ; 


PjriS  :  [  %  1 


Note  that  the  vectors  at  scale  j  +  I  are  half  as  long  as  the  vectors  as  scale  j.  To 
illustrate  further,  suppose  the  wavelet  being  used  is  the  four  coefficient  D4  wavelet, 
and  suppose  one  wants  to  project  from  8  scaling  function  coefficients  at  scale  j  to  4 
scaling  function  coefficients  at  scale  j  -I-  1  and  4  wavelet  coefficients  at  scale  j  -f  1. 


The  decomposition  matrix  in  this  case  is, 

hi  hj  ^3  0  0  0  0 

0  0  hi  h^  0  0 

0  0  0  0  hi  hj  ^3  /l4 

pjJ+i  _  h3  ^4  0  0  0  0  hi  hi  .  . 

~  ^1  53  fll3  5f4  0  0  0  0  ’  ^ 

0  0  5i  52  53  54  0  0 

0  0  0  0  5i  52  53  54 

.  53  54  0  0  0  0  5i  52  . 

where  the  periodicity  is  seen  from  the  coefficients  ‘wrapping  around’. 

Now  let  us  consider  differentiation.  Let  the  four  matrices  ^nxn^  ^NxNi 

and  RjvxJVj  [7]  and  [1],  contain  the  derivative  projection  coefficients, 

:  dj  — »  dj, 

:  Si  I, 

:dj  —*  Sj, 

R^:Sj-^  Sj, 

T 

where  Sj  and  dj  denote  the  coefficients  of  the  expansion  of  the  derivative  of  a  function 
which  is  initially  defined  by  the  expansion  coefficients  Sj  and  dj.  The  exact  form  of  the 
matrices  A,  B,  and  C  is  not  important  for  the  discussion  here.  The  important  point  is 
the  form  of  the  matrix  R.  It  is  always  a  finite  difference  operator.  For  the  D4  wavelet 
R  corresponds  to  the  optimal  central  4th-order  finite  difference  operator.  For  higher 
order  wavelets,  D%,  D%,  etc.,  is  a  finite  difference  operator,  but  it  is  not  optimal 
in  the  sense  of  using  the  minimum  number  of  coefficients  to  obtain  a  given  accuracy. 
The  numerical  values  of  the  coefficients  were  found  in  [1]  and  the  superconvergence 
accuracy  was  proven  in  general  in  [7].  For  the  D4  wavelet  an  explicit  8x8  example 
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of  matrix  R  is, 


One  can  calculate  the  derivative  of  a  wavelet  expansion  at  any  level  in  the  wavelet 
decomposition.  This  subsection  will  explore  the  first  of  three  of  the  options.  To  be 
explicit,  suppose  that  a  periodic  function  f{x)  has  been  approximated  on  a  grid  with 
16  scaling  function  coefficients  to  get  5*0,  and  for  the  current  argument  assume  that 
the  coefficients  have  been  calculated  exactly.  Note  that  periodicity  of  f{x)  induces 
periodicity  on  the  coefficients  Sq.  The  coefficients  of  the  expansion  of  ^/(x)  in  Vo 
are  found  from  sq  by  an  application  of  the  matrix  R^^xie- 
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Let  us  now  examine  the  entire  process  of  going  from  point  values  in  the  physical 
space  to  scaling  function  coefficients  in  Vq,  differentiating,  and  finally  returning  to 
point  values  of  the  differentiated  function  in  the  physical  space.  Begin  with  f{x)  € 
L'^{R)  defined  at  16  evenly-spaced  points  over  [— tt,  v)  and  let  f{x)  is  2?r  periodic. 
To  differentiate  the  samples  of  /(x),  /,  with  the  4-th  order  optimal  central  finite 
difference  operator,  say  we  get, 

l  =  (25) 

Now,  suppose  that  we  have  mapped  these  16  samples  into  the  scaling  function  coeffi¬ 
cients  in  Vq  by  applying  the  circular,  periodicity  implies  circularity,  see  [7],  quadrature 
matrix  Q, 

So  =  Quxiif-  (26) 

^  - 

We  now  find  io  by  applying  ^Aioxie*  The  two  matrices  R  and  Q  are,  however,  both 
circiilar  and,  hence,  commute: 

^0  =  Q\6xl6-^R\6xlBf  >  (27) 

Now,  returning  to  the  physical  space  we  get, 

/  =  QiixieQi^xu-^Riexwfy  (28) 

and  we  are  back  to  equation  (25)  again  since, 

0fd4  =  -^Ri6x16-  (29) 

Hence,  we  have  shown  that  under  the  assumption  of  periodicity  and  without 

data  compression  that  the  D4  wavelet  differentiation  corresponds  exactly  to  optimal 
central  4th-order  finite  differencing.  Note  that  data  compression  is  the  goal  of  any 
wavelet  method.  The  embodiment  of  data  compression  in  the  physical  space  is  a 
nonuniform  grid.  That  is,  the  grid  must  be  dense  in  regions  where  small  structure 
requires  fine  resolution  and  the  grid  can  be  sparse  when  the  data  is  composed  of  large 
scale  components. 


12 


Now  we  move  to  the  first  decomposition  of  Vo  =  VV'j  ©V^i  in  which  data  compression 


can  be  2M:hieved. 

3.3  Wavelet  Expansion  and  Derivative  in  Wi  ©  Vi 

Consider  now  a  decomposition  of  the  vector  of  scaling  function  coefficients  so  onto 
the  scaling  function  and  wavelet  coefficients  at  scale  j  =  1  by  an  application  of  the 
matrix 
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d\ 

.  •*16  - 

L4  J 

(30) 


As  in  Vo,  we  have  16  basis  functions,  but  now  the  subspace  Vo  is  decomposed  into 
‘low  frequency’,  Vi,  and  ‘high  frequency’,  Wi,  components:  Vq  =  Vi  ®  W\.  In  order 
to  calculate  the  coefficients  of  the  derivative  expansion  in  Vi  ©  Wi  the  following 
projections  are  calculated: 

^X8  ^8x8 


*1 

1 

2Aa: 

^8X8  ■^8x8 


*1 

. 

(31) 


If  one  now  applies  the  matrix  (Pfexie)^  {1'  denotes  transpose  and  hence  inverse  for 
this  unitary  matrix)  to  the  derivative  coefficients  at  scale  j  =  1  one  gets, 


[  ^0  ]  —  (^lexie)^ 


Si 


(32) 


and  one  gets  exactly  the  same  coeflBcients  eis  before  when  the  matrix  ^f2?6xi6 
applied  to  sq. 
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was 


Now  suppose  that  f{x)  is  smooth  enough  such  that  a  grid  of  eight  points  provides 
sufficient  resolution.  Define  /a  to  be  the  8  element  vector  containing  every  other  entry 
of  the  16  element  vector  /.  4-th  order  differentiation  of  /a  is  performed  by  applying 

aS^xs  j 


/2  = 


1 


2Ax 


f^xs/a- 


(33) 


Similar  to  above,  we  project  the  eight  dimensional  /a  into  the  eight  dimensional 
wavelet  subspace  Vi  using  Qgxs  and  differentiate  to  get, 


1 

2Ax 


R8x8Qsxsf2y 


(34) 


followed  by  projection  back  into  the  physical  space  with  the  matrix  Qg^g  we  get 
equation  (33)  again. 

That  is,  we  have  seen  that  if  we  work  only  in  Vo  that  we  have  4th-order  finite 
differencing  with  a  grid  spacing  of  Ax,  whereas  if  we  work  only  in  V  we  have  4th- 
order  finite  differencing  with  a  grid  spacing  of  2Ax.  But,  the  two  subspaces  Vo  and  Vi 
are  related  by  Vo  =  Vi  0  Wi .  Recall,  that  the  subspace  Wi  contains  bases  functions 
which  are  locally  oscillatory  and  are  compactly  supported.  An  inner  product  of  these 
bases  with  the  data  /(x)  will  detect  local  oscillations  in  /(x)  and  provide  exactly  the 
information  necessary  to  refine  the  grid  locally  from  2 Ax  to  Ax.  This  wavelet  grid 
refinement  mechanism  can  be  used  not  only  to  add  wavelet  bases  functions  where  one 
has  a  large  inner  product  but  also  to  refine  the  grid  in  the  same  region  and  at  a  scale 
corresponding  to  the  wavelet  scale. 


3.4  Wavelet  Expansion  and  Derivative  in  VTi  0  Vr2  0  ®  Vi 


Let  us  close  this  section  with  a  wavelet  decomposition  that  one  might  use  in  practice. 
That  is,  again  Vo  denotes  the  finest  scale  subspace  and  we  decompose  Vo  as, 


Vo  =  w,eW2eWge  1/3. 


(35) 


14 


The  vector  of  coefficients  in  this  subspace  is  obtained  by  the  following  decompositions: 


Let  us  suppose  that  we  have  performed  the  above  wavelet  decomposition  on  a 
vector  of  data  points,  /,  at  some  point  during  a  simulation  which  contains  data 
at  many  different  scales.  Furthermore,  let  there  be  a  shock,  or  a  near  shock,  near 
the  right-hand  boundary.  The  coefficients  Sj  and  si  represent  local  averages  in  the 
subspace  V3  corresponding  to  the  ‘b2U3e  grid’  of  size  8Ax  and  will  not  yield  much  useful 
information  with  respect  to  the  shock.  The  coefficients  and  d|  of  the  subspace  W3 
will  yield  the  presence  of  oscillations  of  relatively  large  scale.  A  true  shock  contains 
all  frequencies  and  one  would  expect  to  have  some  coefficient  perturbation  even  in 
Wz  yielding  grid  refinement  to  a  grid  spacing  of  4Ai  in  a  neighborhood  of  the  shock. 
The  coefficients  df,  dl,  dg,  and  d|  in  the  subspace  W2  will  detect  oscillations  at  the 
corresponding  scale  only  near  the  shock.  That  is,  the  first  two  coefficients  dj  and  d| 
are  responsible  for  detecting  small  scale  structure  at  the  left-hand  side  of  the  domain 
which  is  away  from  the  shock,  and  we,  therefore,  expect  that  they  will  be  near  zero 
in  magnitude.  The  coefficients  dg  and  d^,  on  the  other  hand,  are  positioned  near 
the  shock  and  will  have  a  relatively  large  amplitude  indicating  the  presence  of  small 
scales.  The  grid  will,  therefore,  be  refined  to  a  spacing  of  2Ax  at  the  right-hand  side 
of  the  domain.  Likewise,  the  remaining  coefficients  in  the  subspace  will  refine  the 
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grid  to  a  spacing  of  Ax  at  the  right-hand  side  of  the  domain. 

In  conclusion,  this  section  has  been  devoted  to  first  illustrating  how  Daubechies- 
based  wavelet  methods  are  in  essence  finite  difference  methods  with  grid  refinement, 
and  second  to  illustrating  how  the  Daubechies-based  wavelets  can  be  used  to  de¬ 
fine  a  grid  for  finite  difference  methods.  The  next  section  will  make  this  symbiotic 
relationship  between  Daubechies  wavelets  and  finite  difference  methods  formal. 
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4  The  Wavelet- Optimized,  Adaptive  Grid  Finite 
Difference  Method,  (WOFD) 

The  new  method  is  to  apply  finite  diflference  on  a  grid  which  is  defined  by  the  magni¬ 
tude  of  wavelet  coefficients  at  various  scales.  That  is,  wavelets  can  detect  oscillations 
in  a  function  at  any  location  and  scale.  Given  a  function  f{x)  for  i  €  /,  where  /  is 
some  interval,  one  decomposes  f(x)  into  a  set  of  wavelet  coefficients  which  depend 
on  two  parameters,  one  for  location  and  one  for  scale,  say  d^,  where  k  is  the  location 
parameter  and  j  is  the  scale  parameter.  If  a  wavelet  coefficient  is  large  in  magnitude, 

141  >  T,  (37) 

or  large  in  energy  (In  practice  the  two  criteria  yield  roughly  the  same  grid.), 

(4)»  >  T,  (38) 

where  T  is  a  coefBcient  threshold  chosen  by  the  user,  then  WOFD  adds  a  grid  point, 
or  two,  at  location  k  and  at  a  grid  density  corresponding  to  the  scale  j.  That  is, 
WOFD  defines  a  grid  which  will  completely  resolve  a  function  across  the  entire  domain 
without  over  resolving  it  where  it  is  relatively  smooth,  or  composed  only  of  large  scale 
structure.  For  the  specific  case  of  the  D4  wavelet  outlined  in  the  previous  section,  the 
D4  wavelet  decomposition  provides  the  optimal  grid  for  4th-order  finite  differencing. 

The  grid  definition  should  be  made  by  a  Daubechies  wavelet  which  corresponds  in 
terms  of  superconvergence  accuracy  to  the  accuracy  of  the  finite  difference  operator. 
That  is,  it  was  proven  in  [7]  that  the  differentiation  matrix  for  the  Daubechies  wavelet 
D2M1  where  M  is  the  number  of  vanishing  moments,  displays  differentiation  accurax:y 
of  order  2M  under  the  assumptions  of  periodicity  and  a  uniform  grid.  Recall,  that  this 
wavelet  subspace  can  only  represent  exsw:tly  the  first  M  polynomials  as  determined 
by  the  number  of  vanishing  moments.  This  order  of  accuracy  2M  should  equal  the 
order  of  accuracy  of  the  finite  difference  operator  for  optimal  grid  selection. 

In  the  next  section  WOFD  will  be  applied  to  Burgers’  equation. 


17 


5  WOFD  Applied  to  Burgers’  Equation 


In  this  section  WOFD  will  be  applied  to  Burgers’  equation, 

dt  dx  ^  dx^  ’ 


(39) 


with  the  initial  condition, 

u(x,  0)  =  i  ^  sin(2;ri).  (40) 

The  goal  of  this  section  is  to  illustrate  that  WOFD  using  the  D4  wavelet  produces 
a  solution  on  a  nonuniform  reduced  grid  which  is  ‘equivalent  in  character’  to  the 
solution  provided  by  4th  order  finite  differencing  on  the  finest  uniform  grid.  That 
is,  for  a  given  viscosity,  e,  there  exists  a  grid  size  fine  enough  such  that  oscillations 
do  not  develop  at  the  ‘shock’.  This  can  be  made  more  precise  by  saying  that  one 
has  a  grid  fine  enough  such  that  the  total  variation  of  the  solution  does  not  increase. 
‘Equivalent  in  character’  means  that  the  total  variation  of  the  solution  provided  by 
WOFD  increases  if  and  only  if  the  total  variation  of  the  solution  produced  by  finite 
differencing  on  the  finest  uniform  grid  increases. 

In  all  the  following  plots  the  uniform  finite  differencing  is  provided  by  the  optimal 
central  4th-order  finite  difference  operator.  The  temporal  discretization  is  achieved 
by  4th-order  Runge-Kutta.  The  WOFD  coefficient  threshold  which  determines  which 
grid  points  to  use  based  on  the  wavelet  coefficient  magnitude  is  set  to  T  =  .001. 
Note  that  when  the  WOFD  coefficient  threshold  is  set  to  T  =  0  that  one  gets  finite 
differencing  on  the  uniform  finest  grid.  In  addition,  if  the  coefficient  threshold  is  set 
to  a  very  large  number,  say  T  =  100,  then  one  gets  finite  differencing  on  a  uniform 
sparse  grid.  The  size  of  this  sparse  grid  is  determined  by  the  number  of  wavelet 
decompositions  one  specifies.  . 


5.1  Periodic  Boundary  Conditions 

In  figure  (1)  WOFD  is  compared  to  finite  differencing  on  uniform  grid  sizes  32,  64, 
and  128.  The  upper  left-hand  plot  has  the  WOFD  solution  superimposed  on  the  finite 
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difference  solution.  The  two  solutions  are  visually  indistinguishable.  Along  the  x-axis 
of  the  plot  an  ‘x’  is  placed  at  every  position  where  a  WOFD  grid  point  is  used.  One 
can  see  that  the  grid  points  are  dense  at  the  shock  and  sparse  where  the  solution  is 
smooth. 

The  remaining  three  plots  show  finite  difference  solutions  on  various  uniform  grid 
sizes.  One  sees  oscillations  for  grid  sizes  32  and  64  but  not  for  grid  size  128. 

Figure  (4)  provides  an  additional  plot  for  periodic  boundary  conditions  with  a 
slightly  larger  viscosity. 

5.2  Nonperiodic  Boundary  Conditions 

The  boundary  conditions  considered  here  are  such  that  the  boundary  values  of  the 
solution  are  required  to  be  fixed  at  the  initial  condition  values, 

u(0,<)  =  u(l,t)  =  (41) 

In  all  the  plots,  differentiation  at  the  boundary  for  the  uniform  finite  difference  method 
is  achieved  by  the  optimal  4th  order  one  sided  finite  difference  coefficients,  see  [2]. 
For  WOFD  both  4th  order  and  3rd  order  boundary  differentiation  will  be  considered. 

In  Figure  (2)  the  differentiation  ai,  the  boundary  is  4th  order,  and,  as  above, 
WOFD  provides  a  solution  which  is  ‘equivalent  in  character’  to  the  finite  difference 
solution  on  the  finest  grid,  128  grid  points,  while  reducing  the  number  of  degrees-of- 
freedom  necessary  to  achieve  this  solution. 

In  Figure  (3)  the  differentiation  at  the  boundary  is  3rd  order.  The  solution  for  the 
3rd-order  boundary  differentiation  is  good,  but  a  slight  difference  can  be  seen  with 
the  finite  difference  method  on  the  finest  grid.  Again,  finite  difference  on  more  coarse 
grids  oscillates  more  at  the  shock  than  the  finest  grid  solution. 

Figure  (5)  provides  an  additional  plot  illustrating  the  solution  provided  by  WOFD 
for  the  nonperiodic  case. 
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FD  on  Grid  128  and  WORD 


FD  on  Grid  128 


Figure  1:  Illustration  of  equivalence  of  WOFD  to  an  equivalent  order  finite  difference 
method  applied  across  the  entire  domain  at  the  finest  scale.  The  boundary  conditions 
are  periodic.  Final  time  =  2,  Viscosity  =  .02, 
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Figure  2:  Illustration  of  equivalence  of  WOFD  to  an  equivalent  order  finite  difference 
method  applied  across  the  entire  domain  at  the  finest  scale.  The  boundary  values 
are  fixed  at  the  initial  condition  values.  Differentiation  at  the  boundary  is  4th  order. 
Final  time  =  .3,  Viscosity  =  .005. 
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FDonG(id128andW0FD 


FD  on  Grid  128 


Figure  3:  Illustration  of  equivalence  of  WOFD  to  an  equivalent  order  finite  difference 
method  applied  across  the  entire  domain  at  the  finest  scale.  The  boundary  values 
are  fixed  at  the  initial  condition  values.  Differentiation  at  the  boundary  is  3rd  order. 
Final  time  =  .3,  Viscosity  =  .005. 
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Figure  5:  WOFD  applied  to  Burgers 
condition  values,  final  time  is  .3,  viscosity 


6  Accuracy  of  WOFD 

In  this  section  the  order  of  accuracy  will  be  examined.  For  numerical  methods  where 
the  grid  is  uniform  the  order  of  accuracy  is  clearly  defined.  For  WOFD,  on  the  other 
hand,  the  discussion  of  accuracy  is  slightly  more  complicated.  That  is,  it  can  be  said 
that  WOFD  approximates  a  4th-order  finite  difference  method  as  well  as  one  desires, 
and  when  the  coefficient  threshold  is  set  to  zero  then  WOFD  is  truly  4th-order.  So, 
WOFD  approximates  methods  of  a  given  order  as  well  as  is  desired.  In  addition,  it 
will  be  seen  that  the  WOFD  method  has  a  very  nice  feature  that  the  rate  of  growth 
of  the  error  in  approximating  the  derivative  is  at  most  a  linear  function  of  frequency. 

6.1  Error  in  Derivative  Approximation 

The  finite  difference  equations  used  in  this  paper  use  either  a  4-point  stencil  or  a 
5-point  stencil,  '"'he  derivative  approximation  error  for  the  equations  on  a  3-point 
stencil  will  be  given.  The  derivative  approximation  error  for  a  larger  stencil  is  an 
obvious  extension  of  the  error  given  here. 

Consider  the  following  Lagrangian  interpolation  of  a  quadratic  polynomial  through 
the  three  points:  (a;i,/(xi)),  (x2>/(®2))>  and  (xa, /(xa)),  for  Xi  <  xa  <  xa: 

Sr(x)  =  (42) 

\  {x-X2){x-X3)  (x-xi)(x-xa)  ,  (x-Xi)(x-X2) 

‘^X,  -  X2)(xi  -  Xa)  ^^^^X2  -  X1)(X2  -  xa)  ^^^\xa  -  x,)(xa  -  xa)' 

If  we  differentiate  ^(x)  and  evaluate  at  xa  we  get, 

(43) 

for  some  a  6  [xj,  xa],  where  Ai  =  xa  —  xi  and  Aa  =  xa  —  xa. 

6.2  Control  of  Error  Growth 

As  given  above  we  will  examine  the  special  case  of  a  3-point  stencil  where  the  grid  is 
refined  by  the  Haar  wavelet.  In  practice  I  never  use  the  Haar  wavelet,  but  it  is  very 
useful  as  an  illustration  tool. 
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As  given  above  the  error  is, 


Err  =  A,A4/"^(a).  (44) 

for  some  a  6  [xi,  X3].  Now,  let  /  be  a  pure  sinusoid  of  frequency  f{x)  =  e‘^*.  Then 
the  magnitude  of  the  error  becomes, 

6|£?rr|  =  AjAj^^  (45) 

The  grid  refinement  mechanism  used  by  WOFD  is  such  that,  roughly, 

A  =  i  (46) 

The  magnitude  of  the  error  becomes, 

6|£rr|  =  (47) 

That  is,  the  refinement  mechanism  keeps  the  rate  of  growth  of  the  error  linear  with 
respect  to  frequency.  Whereas,  without  the  grid  refinement  the  error  grows  as  a  cubic 
in  this  case. 

This  is  one  particular  refinement  mechanism,  but  is  representative  of  a  typical 
refinement  method. 

6.3  Relationship  of  Threshold  Size  to  Solution 

The  grid  for  WOFD  is  chosen  by  the  size  of  wavelet  coefficients  found  from  a  wavelet 
decomposition  of  the  numerical  solution  at  a  given  time.  One  chooses  a  threshold 
with  which  to  measure  the  coefficient  size.  That  is,  if  the  threshold  is  set  to  .001 
then  the  grid  is  refined  at  a  given  location  and  scale  whenever  the  wavelet  coefficients 
at  that  location  and  scale  are  larger  in  magnitude  than  .001.  If  the  threshold  is  set 
to  0  then  one  gets  finite  difference  on  an  evenly-spaced  grid  at  the  finest  scale.  The 
question  then  becomes,  what  is  the  relationship  between  this  threshold  value  and  the 
solution  achieved  by  WOFD.  As  of  now,  a  theoretical  relationship  does  not  exist  but 
a  numerical  relationship  does.  That  is,  if  the  threshold  is  set  to  T  —  le"'’  then  one 


26 


can  expect  that  both  the  li  and  loo  difference  between  WOFD  and  finite  difference 
at  the  finest  scale  will  be  a  constant  times  this  threshold,  say  kT  where  A:  <  10.  For 
example,  for  a  simulation  with  periodic  boundary  conditions,  viscosity  set  at  .01,  and 
the  final  time  set  to  2,  an  l\  difference  of  3.27e~‘‘  and  an  loo  difference  of  2.31e~^ 
were  found  for  a  threshold  value  oi  T  —  le~^.  This  relationship  was  typical  of  all 
simulations  which  were  run. 


7  Stability  of  WOFD 

The  discussion  of  the  stability  of  WOFD  will  be  given  in  terms  of  the  eigenvalues  of 
the  differentiation  matrix. 

7.1  Eigenvalues  of  Differentiation  Matrix 

The  differentiation  matrix  produced  by  this  scheme  will  have  a  full  set  of  eigenvectors. 
We  can,  therefore,  look  at  instability  through  the  magnitude  of  the  eigenvalues. 

Recall  that  the  WOFD  method  can  produce  essentially  a  completely  arbitrary 
grid.  The  differentiation  matrix  can,  therefore,  take  on  an  unlimited  number  of  forms. 
For  this  reason,  an  analytical  approach  is  not  within  reach.  Therefore,  experimental 
results  which  give  the  magnitude  of  the  eigenvalues  of  the  differentiation  after  each 
grid  update  will  be  given.  On  the  following  pages  the  eigenvalues  will  be  given  for 
the  matrix, 

M  =  I  +  VAt  +  1/2D2(a<)2  +  l/6P®(At)3  +  l/24V*{At)\  (48) 

which  corresponds  to  WOFD  being  applied  to  the  linear  equation 

with  4th-order  Runge-Kutta  time  discretization.  The  grid  is  the  grid  that  is  chosen 
for  the  nonlinear  Burgers’  equation.  It  is  seen  that  the  magnitude  of  the  eigenvalues 
do  sometimes  exceed  1,  but  they  rarely  exceed  1  by  very  much.  That  is,  for  the 
periodic  case,  considering  the  maximum  eigenvalue  for  the  4th-order  RK  for  every 
grid  encountered,  the  maximum  eigenvalue  magnitude  for  the  entire  run  up  to  time  2 
is  1.0004.  This  eigenvalue  is  close  enough  to  1  in  magnitude  not  to  excite  instability. 
In  fact,  the  data  would  have  to  have  a  large  component  in  the  direction  of  the  corre¬ 
sponding  eigenvector  and  one  would  have  to  iterate  100  times  to  get  amplification  of 
4%  in  the  direction  of  this  eigenvector. 
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Evais  using  RK4 


Figure  7:  Eigenvalues  of  the  4th  order  Runge-Kutta  differentiation  matrix  at 
The  boundary  values  are  fixed  at  the  initial  condition  values. 
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8  Efficiency  of  WOFD 

In  a  word,  the  efficiency  of  the  WOFD  method  depends  primarily  on  the  data  com¬ 
pression  ratio.  That  is,  we  choose  some  finest  grid  which  captures  all  the  details  of 
our  solution  throughout  the  entire  run.  Let’s  say  that  this  finest  grid  has  N  degrees- 
of-freedom.  Now  choose  how  close,  in  |  •  I2,  you  desire  your  WOFD  solution  to  be 
to  the  solution  on  the  finest  grid.  Choosing  this  ‘closeness  parameter’  dictates  the 
data  compression  ratio.  Let’s  say  that  the  WOFD  method  needs  only  No  degrees- 
of-freedom  to  satisfy  this  criterion.  Then,  the  amount  of  work  done  is,  roughly,  ^ 
times  the  amount  of  work  done  to  get  the  solution  on  the  finest  grid. 

8.1  Work  Involved  for  Grid  Update 

The  grid  update  requires  a  number  of  steps.  I  will  give  a  worst  case  scenario  in 
estimating  the  number  of  operations. 

A  grid  update  requires  order  N  multiplies  where  N  is  the  number  of  degrees-of- 
freedom  in  the  finest  scale.  The  constant  that  is  multiplied  times  N  is  reasonably 
large,  and  the  following  will  show  where  the  operations  are  used: 

1.  One  must  reconstruct  the  function  00  the  finest  grid.  This  requires  about  ION 
multiplies. 

2.  Next,  one  must  perform  a  wavelet  decomposition.  For  a  Daubechies  4  wavelet 
decomposition  the  filters  are  length  4.  Therefore,  the  first  decomposition  re¬ 
quires  AN  multiplies.  Likewise,  the  second  decomposition  requires  2N  multi¬ 
plies.  The  number  of  decompositions  will  determine  the  number  of  multiplies, 
but  let’s  say  that  this  step,  also,  requires  about  ION  multiplies. 

3.  Choosing  a  grid  from  a  wavelet  decomposition  does  not  require  many  operations, 
but  it  does  need  a  number  of  ‘IF-THEN’  statements.  There  is  roughly  1  ‘IF- 
THEN’  statement  for  each  degree-of-freedom.  Let  me,  once  again,  overestimate 
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the  cost  of  each  ‘IF-THEN’  and  also  compensate  for  a  few  other  operations  that 
are  necessary  by  saying  this  requires  roughly  \0N  operations. 

4.  Once  a  grid  is  chosen  the  new  differentiation  matrix  is  found.  This  requires 
about  40iV’o  multiplies  in  a  worst  case  scenario.  The  second  derivative  filter  can 
be  found  from  the  first  derivative  filter  by  a  convolution  of  each  filter.  This 
requires  about  25A/b  multiplies.  Here,  No  is  the  number  grid  points  used  in  the 
compressed  scenario.  Nq  is  some  fraction  of  N. 


The  total  number  of  multiplies  is  the  sum  of  all  the  above  multiples.  That  is, 
about  ZON  +  65No.  These  numbers  are  rough,  and  we  might  as  well  round  up  to 
be  safe  and  say,  50N  +  lOOA/o  multiplies  are  needed  to  define  the  grid  and  build 
a  new  differentiation  matrix.  This  is  reasonably  expensive,  but  the  update  can  be 
done  rarely  during  a  run.  Compare  this  to  finite  difference  on  the  finest  scale.  The 
filters  for  1-st  and  2-nd  order  differentiation  are  length  5.  Each  step  of  Runge-Kutta 
requires,  therefore,  at  least  ION  multiplies.  If  we  are  using  a  fourth  order  RK  then 
we  have  at  least  40iV  mtiltiplies  for  each  time  step.  It  is,  therefore,  fair  to  say  that 
the  grid  update  step  requires  about  the  same  2unount  of  work  as  one  time  step  taken 
using  the  full  grid. 


8.2  Work  Saved  with  Larger  Time  Step 

All  the  numerical  scenarios  use  explicit  time  stepping.  For  Burgers’  equation  with 
viscosity  e  the  time  step  is  set  to, 


At  =  X 


{Axy 


where  Ax  is  the  minimum  spacing  of  the  grid  produce  by  WOFD.  At  the  beginning  of 
any  simulation  if  the  initial  condition  is  smooth,  as  measured  by  a  wavelet  decomposi¬ 
tion,  then  the  minimum  Ax  produced  by  WOFD  is  much  larger  than  the  Ax  used  on 
the  finest  grid.  This  allows  a  much  larger  time  step  without  introducing  large  errors. 
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The  larger  time  step  means  that  far  fewer  total  time  steps  will  need  to  be  taken  to 
arrive  at  the  final  time.  Fewer  time  steps  gives  a  significant  savings  in  terms  of  total 
operations  performed. 
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9  WOFD  in  Higher  Dimensions 


The  effectiveness  of  WOFD  has  been  illustrated  in  1  dimension.  The  natural  follow¬ 
up  question  would  concern  the  effectiveness  of  WOFD  in  higher  dimensions.  At  this 
point  let  us  break  WOFD  into  two  parts:  the  first  part  is  the  grid  definition,  and  the 
second  part  is  the  differencing  on  this  new  grid. 

Grid  definition  falls  within  the  realm  of  local  spectral  analysis.  That  is,  one  is 
interested  in  the  spectrum  locally.  Local  high  frequency  data  requires  a  grid  density 
sufficient  to  resolve  the  highest  frequency,  whereas  local  low  frequency  data  can  be  re¬ 
solved  with  a  relatively  coaurse  grid.  The  wavelet  structure  provides  a  very  convenient 
mechanism  to  perform  this  nested  group  of  short-time  Fourier  transforms.  For  higher 
dimensions,  say  2  dimensions,  one  need  only  choose  a  coordinate  system  for  the  space 
and  simply  perform  the  wavelet  filtering  throughout  all  of  the  data  and  parallel  to 
each  axis.  The  most  common  question  at  this  point  concerns  the  effectiveness  of  this 
method  of  grid  selection  when  the  data  is  composed  of  structure  which  is  at  a  45 
degree  angle  to  the  axes.  The  simple  answer  is  that  it  works  well  since  all  structures 
within  the  data  can  be  projected  onto  the  orthogonal  coordinate  system  which  spans 
the  space.  If,  however,  one  is  not  satisfied  with  the  grid  given  in  this  situation  then 
the  parameter  which  adjusts  the  sensitivity  of  grid  selection  can  be  adjusted.  With 
this  sensitivity  adjustment,  one  will  always  find  a  suitable  grid.  Included  here  are 
three  sets  of  data  which  will  illustrate  the  effective  grid  selection.  The  first  function 
is  a  discontinuity  at  a  45  degree  angle  to  the  axes,  see  Figure  (8),  and  the  grid  chosen 
for  this  data,  see  (9).  The  second  set  of  data  corresponds,  roughly,  to  the  inner  and 
outer  flow  near  a  boundary, 

/(a:,y)  =  1 -e“^,  (49) 

see  Figure  (10),  and  the  grid  chosen  in  this  case,  see  Figure  (11).  The  third  set  of  data 
is  numerically-generated  pressure  from  a  turbulent  flow.  A  contour  diagram  is  given 
in  Figure  (12)  and  the  grid  chosen  is  given  in  Figure  (13).  In  all  cases  the  wavelet  is 
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the  D4  wavelet. 

Differencing  on  a  grid  chosen  by  WOFD  will  depend  on  the  application.  The  D4  is 
the  optimal  wavelet  if  one  is  using  a  central  optimal  4th-order  finite  difference  method 
in  the  sense  of  Section  3.  But  one  is  not  confined  to  matching  the  wavelet  precisely 
to  the  differencing  method.  The  only  recommendation  is  that  the  stencil  of  the  finite 
difference  method  be  of  roughly  the  same  size  as  the  length  of  the  wavelet  filters. 
This  will  insure  that  in  the  grid  selection  one  is  not  filtering  data  which  is  too  far 
outside  of  the  support  of  the  differencing  stencil.  Also,  depending  on  the  application 
one  might  need  a  reliable  mechanism  for  choosing  locally  rectangular  grids.  In  this 
case  one  would  choose  the  finest  grid  produced  by  WOFD  in  each  rectangular  region. 

One  other  issue  is  the  use  of  WOFD  to  refine  beyond  the  ‘finest  grid’.  For  the  1 
dimension  case  considered  in  this  oa'^er  there  was  always  an  underlying  finest  grid. 
This  was  a  theoretical  convenience  and  not  a  necessity  or  even  a  desirable  feature. 
Based  on  the  energy  or  magnitude  of  the  smallest  scale  wavelet  coefficients  one  can 
refine  to  any  desirable  grid  density,  see  [10]. 
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Figure  9:  The  grid  chosen  by  the  wavelet  decomposition  for  a  discontinuity  at  a  45 
degree  angle  to  the  axes. 
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Figure  11:  The  grid  chosen  by  the  wavelet  decomposition  of  the  flow  near  a  boundary. 


10  Conclusion 


The  WOFD  method  is  essentially  the  same  as  a  Daubechies-based  wavelet  method. 
However,  two  serious  problems  encoimtered  with  wavelet  numerical  methods  are  over¬ 
come.  That  is,  boundary  conditions  can  now  be  imposed  in  exactly  the  same  manner 
they  are  imposed  for  finite  difference  methods.  Furthermore,  there  is  no  longer  a 
problem  with  nonlinear  terms  since  we  are  now  working  with  the  point  values  of  the 
function  in  the  physical  space.  WOFD  will  approximate  a  conservative  numerical 
method  as  long  as  one  is  working  with  a  conservative  numerical  method  on  the  finest 
uniform  grid.  This  approximation  can  be  as  fine  as  the  user  wishes. 

In  this  paper  the  WOFD  method  has  been  explored  for  the  case  of  a  5-point  4- 
th  order  finite  difference  operator  on  an  arbitrary  grid.  We  can,  however,  use  this 
method  with  higher  order  schemes  by  using  the  filters  associated  with  the  higher 
order  Daubechies  wavelets  to  define  the  grid.  The  only  suggestion  is  that  the  stencil 
size  of  the  numerical  scheme  be  of  roughly  the  same  size  as  the  length  of  the  filters. 

Higher  dimensional  applications  of  WOFD  remain  to  be  explored.  As  mentioned, 
WOFD  can  be  broken  into  two  parts,  the  grid  selection  and  the  differencing.  The 
grid  selection  step  requires  the  wavelet  analysis  to  detect  the  local  oscillation  content 
of  the  data.  Based  on  this  grid,  one  can  choose  a  number  of  ways  to  apply  finite 
differencing.  It  has  been  shown  here  that  grid  selection  is  quite  effective  for  ‘difficult’ 
data  sets  in  two  dimensions.  Future  research  will  combine  this  2-dimensionaI  grid 
selection  with  appropriate  differencing. 
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